Learning Scalable Deep Kernels with Recurrent Structure

Many applications in speech, robotics, finance, and biology deal with sequential data, where ordering matters and recurrent structures are common. However, this structure cannot be easily captured by standard kernel functions. To model such structure, we propose expressive closed-form kernel functions for Gaussian processes. The resulting model, GP-LSTM, fully encapsulates the inductive biases of long short-term memory (LSTM) recurrent networks, while retaining the non-parametric probabilistic advantages of Gaussian processes. We learn the properties of the proposed kernels by optimizing the Gaussian process marginal likelihood using a new provably convergent semi-stochastic gradient procedure, and exploit the structure of these kernels for scalable training and prediction. This approach provides a practical representation for Bayesian LSTMs. We demonstrate state-of-the-art performance on several benchmarks, and thoroughly investigate a consequential autonomous driving application, where the predictive uncertainties provided by GP-LSTM are uniquely valuable.


Introduction
There exists a vast array of machine learning applications where the underlying datasets are sequential.Applications range from the entirety of robotics, to speech, audio and video processing.While neural network based approaches have dealt with the issue of representation learning for sequential data, the important question of modeling and propagating uncertainty across time has rarely been addressed by these models.For a robotics application such as a self-driving car, however, it is not just desirable, but essential to have complete predictive densities for variables of interest.When trying to stay in lane and keep a safe following distance from the vehicle front, knowing the uncertainty associated with lanes and lead vehicles is as important as the point estimates.
Recurrent models with long short-term memory (LSTM) (Hochreiter and Schmidhuber, 1997) have recently emerged as the leading approach to modeling sequential structure.The LSTM is an efficient gradient-based method for training recurrent networks.LSTMs use a memory cell inside each hidden unit and a special gating mechanism that stabilizes the flow of the back-propagated errors, improving the learning process of the model.While the LSTM provides state-of-the-art results on speech and text data (Graves et al., 2013;Sutskever et al., 2014), quantifying uncertainty or extracting full predictive distributions from deep models is still an area of active research (Gal and Ghahramani, 2016a).
In this paper, we quantify the predictive uncertainty of deep models by following a Bayesian nonparametric approach.In particular, we propose kernel functions which fully encapsulate the structural properties of LSTMs, for use with Gaussian processes.The resulting model enables Gaussian processes to achieve state-of-the-art performance on sequential regression tasks, while also allowing for a principled representation of uncertainty and non-parametric flexibility.Further, we develop a provably convergent semi-stochastic optimization algorithm that allows mini-batch updates of the recurrent kernels.We empirically demonstrate that this semi-stochastic approach significantly improves upon the standard non-stochastic first-order methods in runtime and in the quality of the converged solution.
For additional scalability, we exploit the algebraic structure of these kernels, decomposing the relevant covariance matrices into Kronecker products of circulant matrices, for O(n) training time and O(1) test predictions (Wilson et al., 2015;Wilson and Nickisch, 2015).Our model not only can be interpreted as a Gaussian process with a recurrent kernel, but also as a deep recurrent network with probabilistic outputs, infinitely many hidden units, and a utility function robust to overfitting.
Throughout this paper, we assume basic familiarity with Gaussian processes (GPs).We provide a brief introduction to GPs in the background section; for a comprehensive reference, see, e.g., Rasmussen and Williams (2006).In the following sections, we formalize the problem of learning from sequential data, provide background on recurrent networks and the LSTM, and present an extensive empirical evaluation of our model.Specifically, we apply our model to a number of tasks, including system identification, energy forecasting, and self-driving car applications.Quantitatively, the model is assessed on the data ranging in size from hundreds of points to almost a million with various signal-to-noise ratios demonstrating state-of-the-art performance and linear scaling of our approach.Qualitatively, the model is tested on consequential self-driving applications: lane estimation and lead vehicle position prediction.Indeed, the main focus of this paper is on achieving stateof-the-art performance on consequential applications involving sequential data, following straightforward and scalable approaches to building highly flexible Gaussian process.
We release our code as a library at: http://github.com/alshedivat/keras-gp.This library implements the ideas in this paper as well as deep kernel learning (Wilson et al., 2016a) via a Gaussian process layer that can be added to arbitrary deep architectures and deep learning frameworks, following the Keras API specification.More tutorials and resources can be found at https://people.orie.cornell.edu/andrew/code.

Background
We consider the problem of learning a regression function that maps sequences to real-valued target vectors.Formally, let X = {x i } n i=1 be a collection of sequences, , each with corresponding length, l i , where x j i ∈ X , and X is an arbitrary domain.Let y = {y i } n i=1 , y i ∈ R d , be a collection of the corresponding real-valued target vectors.Assuming that only the most recent L steps of a sequence are predictive of the targets, the goal is to learn a function, f : X L → R d , from some family, F, based on the available data.
As a working example, consider the problem of estimating position of the lead vehicle at the next time step from LIDAR, GPS, and gyroscopic measurements of a self-driving car available for a number of previous steps.This task is a classical instance of the sequenceto-reals regression, where a temporal sequence of measurements is regressed to the future position estimates.In our notation, the sequences of inputs are vectors of measurements, are indexed by time and would be of growing lengths.Typically, input sequences are considered up to a finite-time horizon, L, that is assumed to be predictive for the future targets of interest.The targets, y 1 , y 2 , . . ., y n , are two-dimensional vectors that encode positions of the lead vehicle in the ego-centric coordinate system of the self-driving car.
Note that the problem of learning a mapping, f : X L → R d , is challenging.While considering whole sequences of observations as input features is necessary for capturing long-term temporal correlations, it virtually blows up the dimensionality of the problem.If we assume that each measurement is p-dimensional, i.e., X ⊆ R p , and consider L previous steps as distinct features, the regression problem will become (L × p)-dimensional.Therefore, to avoid overfitting and be able to extract meaningful signal from a finite amount of data, it is crucial to exploit the sequential nature of observations.
Recurrent models.One of the most successful ways to exploit sequential structure of the data is by using a class of recurrent models.In the sequence-to-reals regression scenario, such a model expresses the mapping f : X L → R d in the following general recurrent form: where x t is an input observation at time t, h t is a corresponding latent representation, and y is a target vector.Functions φ(•) and ψ(•) specify model transitions and emissions, respectively, and δ t and t are additive noises.While φ(•) and ψ(•) can be arbitrary, they are typically time-invariant.This strong but realistic assumption incorporated into the structure of the recurrent mapping significantly reduces the complexity of the family of functions, F, regularizes the problem, and helps to avoid severe overfitting.
Recurrent models can account for various patterns in sequences by memorizing internal representations of their dynamics via adjusting φ and ψ.Recurrent neural networks (RNNs) model recurrent processes by using linear parametric maps followed by nonlinear activations: where W hy , W hh , W xh are weight matrices to be learned 1 and φ(•) and ψ(•) here are some fixed element-wise functions.Importantly and contrary to the standard hidden Markov models (HMMs), the state of an RNN at any time t is distributed and effectively represented by an entire hidden sequence, A major disadvantage of the vanilla RNNs is that their training is nontrivial due to the so-called vanishing gradient problem (Bengio et al., 1994): the error back-propagated through t time steps diminishes exponentially which makes learning long-term relationships nearly impossible.LSTM.To overcome vanishing gradients, Hochreiter and Schmidhuber (1997) proposed a long short-term memory (LSTM) mechanism that places a memory cell into each hidden unit and uses differentiable gating variables.The update rules for the hidden representation at time t have the following form (here σ (•) and tanh (•) are element-wise sigmoid and hyperbolic tangent functions, respectively): (3) As illustrated above, g t i , g t c , and g t o correspond to the input, forget, and output gates, respectively.These variables take their values in [0, 1] and when combined with the internal states, c t , and inputs, x t , in a multiplicative fashion, they play the role of soft gating.The gating mechanism not only improves the flow of errors through time, but also, allows the the network to decide whether to keep, erase, or overwrite certain memorized information based on the forward flow of inputs and the backward flow of errors.This mechanism adds stability to the network's memory.
Gaussian processes.The Gaussian process (GP) is a Bayesian nonparametric model that generalizes the Gaussian distributions to functions.We say that a random function f is drawn from a GP with a mean function µ and a covariance kernel k, f ∼ GP (µ, k), if for any vector of inputs, [x 1 , x 2 , . . ., x n ], the corresponding vector of function values is Gaussian: with mean µ, such that µ i = µ(x i ), and covariance matrix K X,X that satisfies (K X,X ) ij = k(x i , x j ).GPs can be seen as distributions over the reproducing kernel Hilbert space (RKHS) of functions which is uniquely defined by the kernel function, k (Schölkopf and Smola, 2002).
GPs with RBF kernels are known to be universal approximators with prior support to within an arbitrarily small epsilon band of any continuous function (Micchelli et al., 2006).Assuming additive Gaussian noise, y | x ∼ N f (x), σ 2 , and a GP prior on f (x), given training inputs x and training targets y, the predictive distribution of the GP evaluated at an arbitrary test point x * is: where (5) Here, K X * ,X , K X,X * , K X,X , and K X * ,X * are matrices that consist of the covariance function, k, evaluated at the corresponding points, x ∈ X and x * ∈ X * , and µ X * is the mean function evaluated at x * ∈ X * .GPs are fit to the data by optimizing the evidence-the marginal probability of the data given the model-with respect to kernel hyperparameters.The evidence has the form: where we use a shorthand K for K X,X , and K implicitly depends on the kernel hyperparameters.This objective function consists of a model fit and a complexity penalty term that results in an automatic Occam's razor for realizable functions (Rasmussen and Ghahramani, 2001).By optimizing the evidence with respect to the kernel hyperparameters, we effectively learn the the structure of the space of functional relationships between the inputs and the targets.For further details on Gaussian processes and relevant literature we refer interested readers to the classical book by Rasmussen and Williams (2006).
Turning back to the problem of learning from sequential data, it seems natural to apply the powerful GP machinery to modeling complicated relationships.However, GPs are limited to learning only pairwise correlations between the inputs and are unable to account for long-term dependencies, often dismissing complex temporal structures.Combining GPs with recurrent models has potential to addresses this issue.

Related work
The problem of learning from sequential data, especially from temporal sequences, is well known in the control and dynamical systems literature.Stochastic temporal processes are usually described either with generative autoregressive models (AM) or with statespace models (SSM) (Van Overschee and De Moor, 2012).The former approach includes nonlinear auto-regressive models with exogenous inputs (NARX) that are constructed by using, e.g., neural networks (Lin et al., 1996) or Gaussian processes (Kocijan et al., 2005).The latter approach additionally introduces unobservable variables, the state, and constructs autoregressive dynamics in the latent space.This construction allows to represent and propagate uncertainty through time by explicitly modeling the signal (via the state evolution) and the noise.Generative SSMs can be also used in conjunction with discriminative models via the Fisher kernel (Jaakkola and Haussler, 1999).
Modeling time series with GPs is equivalent to using linear-Gaussian autoregressive or SSM models (Box et al., 1994).Learning and inference are efficient in such models, but they are not designed to capture long-term dependencies or correlations beyond pairwise.Wang et al. (2005) introduced GP-based state-space models (GP-SSM) that use GPs for transition and/or observation functions.These models appear to be more general and flexible as they account for uncertainty in the state dynamics, though require complicated approximate training and inference, which are hard to scale (Turner et al., 2010;Frigola et al., 2014).
Perhaps the most recent relevant work to our approach is recurrent Gaussian processes (RGP) (Mattos et al., 2015).RGP extends the GP-SSM framework to regression on sequences by using a recurrent architecture with GP-based activation functions.The structure of the RGP model mimics the standard RNN, where every parametric layer is substituted with a Gaussian process.This procedure allows one to propagate uncertainty throughout the network for an additional cost.Inference is intractable in RGP, and efficient training requires a sophisticated approximation procedure, the so-called recurrent variational Bayes.In addition, the authors have to turn to RNN-based approximation of the variational mean functions to battle the growth of the number of variational parameters with the size of data.While technically promising, RGP seems problematic from the application perspective, especially in its implementation and scalability aspects.
Our model has several distinctions with prior work aiming to regress sequences to reals.Firstly, one of our goals is to keep the model as simple as possible while being able to represent and quantify predictive uncertainty.We maintain an analytical objective function and refrain from complicated and difficult-to-diagnose inference schemes.This simplicity is achieved by giving up the idea of propagating signal through a chain GPs connected in a recurrent fashion.Instead, we propose to directly learn kernels with recurrent structure via joint optimization of a simple functional composition of a standard GP with a recurrent model (e.g., LSTM), as described in detail in the following section.Similar approaches have recently been explored and proved to be fruitful for non-recurrent deep networks (Wilson et al., 2016a,b;Calandra et al., 2016).We remark that combinations of GPs with nonlinear functions have also been considered in the past in a slightly different setting of warped regression targets (Snelson et al., 2003;Wilson and Ghahramani, 2010;Lázaro-Gredilla, 2012).Additionally, uncertainty over the recurrent parts of our model is represented via dropout, which is computationally cheap and turns out to be equivalent to approximate Bayesian inference in a deep Gaussian process (Damianou and Lawrence, 2013) with particular intermediate kernels (Gal and Ghahramani, 2016a,b).Finally, one can also view our model as a standalone flexible Gaussian process, which leverages learning techniques that scale to massive datasets (Wilson and Nickisch, 2015;Wilson et al., 2015).

Learning recurrent kernels
Gaussian processes with different kernel functions correspond to different structured probabilistic models.For example, some special cases of the Matérn class of kernels induce models with Markovian structure (Stein, 1999).To construct deep kernels with recurrent structure we transform the original input space with an LSTM network and build a kernel directly in the transformed space, as shown in Figure 1b.
In particular, let L φ : X L → H be an arbitrary deterministic transformation of the input sequences into some latent space, H. Next, let k : H 2 → R be a real-valued kernel defined on H.The decomposition of the kernel and the transformation, k = k • φ, is defined as It is trivial to show that k(x, x ) is a valid kernel defined on X L (MacKay, 1998, Ch. 5.4.3).
In addition, if φ(•) is represented by a neural network, the resulting model can be viewed as the same network, but with an additional GP-layer and the negative log marginal likelihood (NLML) objective function used instead of the standard mean squared error (MSE).
Input embedding is well-known in the Gaussian process literature (e.g.MacKay, 1998;Hinton and Salakhutdinov, 2008).Recently, Wilson et al. (2016a,b) have successfully scaled the approach and demonstrated strong results in regression and classification tasks for kernels based on feedforward and convolutional architectures.In this paper, we apply the same technique to learn kernels with recurrent structure by transforming input sequences with a recurrent neural network that acts as φ(•).In particular, a (multi-layer) LSTM architecture is used to embed L steps of the input time series into a single vector in the hidden space, H.For the embedding, as common, we use the last hidden vector produced by the recurrent network.Note that however, any variations to the embedding (e.g., using other types of recurrent units, adding 1-dimensional pooling layers, or attention mechanisms) are all fairly straightforward2 .More generally, the recurrent transformation can be random itself (Figure 1), which would enable direct modeling of uncertainty within the recurrent dynamics, but would also require inference for φ (e.g., as in Mattos et al., 2015).In this study, we limit our consideration of random recurrent maps to only those induced by dropout.
Unfortunately, once the the MSE objective is substituted with NLML, it no longer factorizes over the data.This prevents us from using the well-established stochastic optimization techniques for training our recurrent model.In the case of feedforward and convolutional networks, Wilson et al. (2016a) proposed to pre-train the input transformations and then fine-tune them by jointly optimizing the GP marginal likelihood with respect to hyperparameters and the network weights using full-batch algorithms.When the transformation is recurrent, stochastic updates play a key role.Therefore, we propose a semi-stochastic block-gradient optimization procedure which allows mini-batching weight updates and fully joint training of the model from scratch.

Optimization
The negative log marginal likelihood of the Gaussian process has the following form: where and implicitly depends on the base kernel hyperparameters, θ, and the parameters of the recurrent neural transformation, φ(•), denoted W and further referred as the transformation hyperparameters.Our goal is to optimize L with respect to both θ and W .
The derivative of the NLML objective with respect to θ is standard and takes the following form (Rasmussen and Williams, 2006): where ∂K/∂θ is depends on the kernel function, k(•, •), and usually has an analytic form.
The derivative with respect to the l-th transformation hyperparameter, W l , is as follows: ) where h i = φ(x i ) corresponds to the latent representation of the the i-th data instance.Once the derivatives are computed, the model can be trained with any first-order or quasi-Newton optimization routine.However, application of the stochastic gradient method-the de facto standard optimization routine for deep recurrent networks-is not straightforward: neither the objective, nor its derivatives factorize over the data4 due to the kernel matrix inverses, and hence convergence is not guaranteed.
Semi-stochastic alternating gradient descent.Observe that once the kernel matrix, K, is fixed, the expression, K −1 yy K −1 − K −1 , can be precomputed on the full data and fixed.Subsequently, Eq. ( 10) turns into a weighted sum of independent functions of each data point.This observation suggests that, given a fixed kernel matrix, one could compute a stochastic update for W on a mini-batch of training points by only using the corresponding sub-matrix of K. Hence, we propose to optimize GPs with recurrent kernels in a semi-stochastic fashion, alternating between updating the kernel hyperparameters, θ, on the full data first, and then updating the weights of the recurrent network, W , using stochastic steps.The procedure is given in Algorithm 1.
Semi-stochastic alternating gradient descent is a special case of block-stochastic gradient iteration (Xu and Yin, 2015).While the latter splits the variables into arbitrary blocks and applies Gauss-Seidel type stochastic gradient updates to each of them, our procedure alternates between applying deterministic updates to θ and stochastic updates to W of the form θ 5 .The corresponding Algorithm 1 is provably convergent for convex and non-convex problems under certain conditions.The following theorem adapts results of Xu and Yin (2015) to our optimization scheme.
Theorem 1 (informal) Semi-stochastic alternating gradient descent converges to a fixed point when the learning rate, λ t , decays as Θ(1/t 1+δ 2 ) for any δ ∈ (0, 1]. Applying alternating gradient to our case has a catch: the kernel matrix (and its inverse) has to be updated each time W and θ are changed, i.e., on every mini-batch iteration (marked red in Algorithm 1).Computationally, this updating strategy defeats the purpose of stochastic gradients because we have to use the entire data on each step.To deal with the issue of computational efficiency, we use ideas from asynchronous optimization.

4:
for all mini-batches X b in X do 5: w ← w + update w (X b , w, K stale ).Update the kernel matrix, K. 8: until Convergence output Optimal θ * and w * Asynchronous techniques.One of the recent trends in parallel and distributed optimization is applying updates in an asynchronous fashion (Agarwal and Duchi, 2011).Such strategies naturally require some tolerance to delays in parameter updates (Langford et al., 2009).In our case, we modify Algorithm 1 to allow delayed kernel matrix updates.
The key observation is very intuitive: when the stochastic updates of W are small enough, K does not change much between mini-batches, and hence we can perform multiple stochastic steps for W before re-computing the kernel matrix, K, and still converge.For example, K may be updated once at the end of each pass through the entire data (see Algorithm 2).To ensure convergence of the algorithm, it is important to strike the balance between (a) the learning rate for W and (b) the frequency of the kernel matrix updates.
The following theorem provides convergence results under certain conditions.
Formal statements, conditions, and proofs for Theorems 1 and 2 are given in Appendix C.2.
Why stochastic optimization?GPs with recurrent kernels can be also trained with full-batch gradient descent, as proposed by Wilson et al. (2016a).However, stochastic gradient methods have been proved to attain better generalization (Hardt et al., 2016) and often demonstrate superior performance in deep and recurrent architectures (Wilson and Martinez, 2003).Moreover, stochastic methods are 'online', i.e., they update model parameters based on subsets of an incoming data stream, and hence can scale to very large datasets.In our experiments, we demonstrate that GPs with recurrent kernels trained with Algorithm 2 converge faster (i.e., require fewer passes through the data) and attain better performance than if trained with full-batch techniques.
Stochastic variational inference.Stochastic variational inference (SVI) in Gaussian processes (Hensman et al., 2013) is another viable approach to enabling stochastic optimization for GPs with recurrent kernels.Such method would optimize a variational lower bound on the original objective that factorizes over the data by construction.Recently, Wilson et al. (2016b) developed such a stochastic variational approach in the context of deep kernel learning.Note that unlike all previous existing work, our proposed approach does not require a variational approximation to the marginal likelihood to perform mini-batch training of Gaussian processes.

Scalability
Learning and inference with Gaussian processes requires solving a linear system involving an n × n kernel matrix, K −1 y, and computing a log determinant over K.These operations typically require O(n 3 ) computations for n training data points, and O(n 2 ) storage.In our approach, scalability is achieved through semi-stochastic training and structure-exploiting inference.In particular, asynchronous semi-stochastic gradient descent reduces both the total number of passes through the data required for the model to converge and the number of calls to the linear system solver; exploiting the structure of the kernels significantly reduces the time and memory complexities of the linear algebraic operations.
More precisely, we replace all instances of the covariance matrix K y with W K U,U W , where W is a sparse interpolation matrix, and K U,U is the covariance matrix evaluated over m latent inducing points, which decomposes into a Kronecker product of circulant matrices (Wilson and Nickisch, 2015;Wilson et al., 2015).This construction makes inference and learning scale as O(n) and test predictions be O(1), while preserving model structure.For the sake of completeness, we provide an overview of the underlying algebraic machinery in Appendix A.
At a high level, because W is sparse and K U,U is structured it is possible to take extremely fast matrix vector multiplications (MVMs) with the approximate covariance matrix K X,X .One can then use methods such as linear conjugate gradients, which only use MVMs, to efficiently solve linear systems.MVM or scaled eigenvalue approaches (Wilson and Nickisch, 2015;Wilson et al., 2015) can also be used to efficiently compute the log determinant and its derivatives.Kernel interpolation (Wilson et al., 2015) also enables fast predictions, as we describe further in the Appendix.

Experiments
We compare the proposed Gaussian processes with recurrent kernels based on RNN and LSTM architectures (GP-RNN/LSTM) with a number of baselines on datasets of various complexity and ranging in size from hundreds to almost a million of time points.For the datasets with more than a few thousand points, we use a massively scalable version of GPs (see Section 4.2) and demonstrate its scalability during inference and learning.We carry

Data and the setup
Below, we describe each dataset we used in our experiments and the associated prediction tasks.The essential statistics for the datasets are summarized in Table 1.System identification.In the first set of experiments, we used publicly available nonlinear system identification datasets: Actuator 6 (Sjöberg et al., 1995) and Drives 7 (Wigren, 2010).Both datasets had one dimensional input and output time series.Actuator had the size of the valve opening as the input and the resulting change in oil pressure as the output.Drives was from a system with motors that drive a pulley using a flexible belt; the input was the sum of voltages applied to the motors and the output was the speed of the belt.
Smart grid data 8 .We considered the problem of forecasting for the smart grid that consisted of two tasks (Figure 2).The first task was to predict power load from the historical temperature data.The data had 11 input time series coming from hourly measurements of temperature on 11 zones and an output time series that represented the cumulative hourly power load on a U.S. utility.The second task was to predict power generated by wind farms from the wind forecasts.The data consisted of 4 different hourly forecasts of the wind and hourly values of the generated power by a wind farm.Each wind forecast was a 6. http://www.iau.dtu.dk/nnbook/systems.html7. http://www.it.uu.se/research/publications/reports/2010-020/NonlinearData.zip.8.The smart grid data were taken from Global Energy Forecasting Kaggle competitions organized in 2012.4-element vector that corresponded to zonal component, meridional component, wind speed and wind angle.In our experiments, we concatenated the 4 different 4-element forecasts, which resulted in a 16-dimensional input time series.
Self-driving car dataset9 .One of the main target applications of the proposed model is prediction for autonomous driving.We considered a large dataset coming from sensors of a self-driving car that was recorded on two car trips with discretization of 10 ms.The data featured two sets of GPS ECEF locations, ECEF velocities, measurements from a fiber-optic gyro compass, LIDAR, and a few more time series from a variety of IMU sensors.Additionally, locations of the left and right lanes were extracted from a video stream for each time step as well as the position of the lead vehicle from the LIDAR measurements.We considered the data from the first trip for training and from the second trip for validation and testing.A visualization of the car routes with 25 second discretization in the ENU coordinates are given in Figure 3.We consider four tasks, the first two of which are more of proof-of-concept type variety, while the final two are fundamental to good performance for a self-driving car: 1. Speed prediction from noisy GPS velocity estimates and gyroscopic inputs.
2. Prediction of the angular acceleration of the car from the estimates of its speed and steering angle.
3. Point-wise prediction of the lanes from the estimates at the previous time steps, and estimates of speed, gyroscopic and compass measurements.
4. Prediction of the lead vehicle location from its location at the previous time steps, and estimates of speed, gyroscopic and compass measurements.We provide more specific details on the smart grid data and self-driving data in Appendix D. Models and metrics.We used a number of classical baselines: NARX (Lin et al., 1996), GP-NARX models (Kocijan et al., 2005), and classical RNN and LSTM architectures.The kernels of our models, GP-NARX/RNN/LSTM, used the ARD base kernel function and were structured by the corresponding baselines. 10As the primary metric, we used root mean squared error (RMSE) on a held out set and additionally negative log marginal likelihood (NLML) on the training set for the GP-based models.
We train all the models to perform one-step-ahead prediction in an autoregressive setting, where targets at the future time steps are predicted from the input and target values at a fixed number of past time steps.For the system identification task, we additionally consider the non-autoregressive scenario (i.e., mapping only input sequences to the future targets), where we are performing prediction in the free simulation mode, and included recurrent Gaussian processes (Mattos et al., 2015) in our comparison.In this case, none of the future targets are available and the models have to re-use their own past predictions to produce future forecasts).
A note on implementation.Recurrent parts of each model were implemented using Keras 11 library.We extended Keras with the Gaussian process layer and developed a backed engine based on the GPML library 12 .Our approach allows us to take full advantage of the functionality available in Keras and GPML, e.g., use automatic differentiation for the recurrent part of the model.Our code is available at http://github.com/alshedivat/kgp/.The inset zooms the region of the plot right beneath it using log scale for the vertical axis.full and mini denote full-batch and mini-batch optimization procedures, respectively, while 16 and 64 refer to models with the respective number of units per hidden layer.Two charts on the right: Test RMSE for a given architecture trained with a specified method and/or learning rate.

Analysis
This section discusses quantitative and qualitative experimental results.We only briefly introduce the model architectures and the training schemes used in each of the experiments.We provide a comprehensive summary of these details in Appendix E.

Convergence of the optimization
To address the question of whether stochastic optimization of recurrent kernels is necessary and to assess the behavior of the proposed optimization scheme with delayed kernel updates, we conducted a number of experiments on the Actuator dataset (Figure 4).First, we constructed two GP-LSTM models with 1 recurrent hidden layer and 16 or 64 hidden units and trained them with (non-stochastic) full-batch iterative procedure (similar to the proposal of Wilson et al. (2016a)) and with our semi-stochastic optimizer with delayed kernel updates (Algorithm 2).The convergence results are given on the first two charts.Both in terms of the error on a held out set and the NLML on the training set, the models trained with mini-batches converged faster and demonstrated better final performance.
Next, we compared the two optimization schemes on the same GP-LSTM architecture with different sizes of the hidden layer ranging from 2 to 32.It is clear from the third chart that, even though full-batch approach seemed to find a better optimum when the number of hidden units was small, the stochastic approach was clearly superior for larger hidden layers.
Finally, we compared the behavior of Algorithm 2 with different number of mini-batches used for each epoch (equivalently, the number of steps between the kernel matrix updates) and different learning rates.The results are give on the last chart.As expected, there is a fine balance between the number of mini-batches and the learning rate: if the number of mini-batches is large (i.e., the delay between the kernel updates becomes too long) while the learning rate is high enough, optimization does not converge; at the same time, an appropriate combination of the learning rate and the mini-batch size leads better generalization than the default batch approach of Wilson et al. (2016a).

Regression, Auto-regression, and Free Simulation
In this set of experiments, our main goal is to provide a comparison between three different modes of one-step-ahead prediction, referred to as (i) regression, (ii) autoregression, and (iii) free simulation, and compare performance of our models with RGP-a classical RNN with every parametric layer substituted with a Gaussian process (Mattos et al., 2015)-on the Actuator and Drives datasets.The difference between the prediction modes consists in whether and how the information about the past targets is used.In the regression scenario, inputs and targets are separate time series and the model learns to map input values at a number of past time points to a target value at a future point in time.Autoregression, additionally, uses the true past target values as inputs; in the free simulation mode, the model learns to map past inputs and its own past predictions to a future target.
In the experiments in autoregression and free simulation modes, we used short time lags, L = 10, as suggested by Mattos et al. (2015).In the regression mode, since the model does not build the recurrent relationships based on the information about the targets (or their estimates), it generally requires larger time lags that can capture the state of the dynamics.Hence we increased the time lag to 32 in the regression mode.More details are given in Appendix E.
We present the results in Table 3.We note that GP-based architectures consistently yielded improved predictive performance compared to their vanilla deep learning counterparts on both of the datasets, in each mode.Given the small size of the datasets, we attribute such behavior to better regularization properties of the negative log marginal likelihood loss function.We also found out that when GP-based models were initialized with weights of pre-trained neural networks, they tended to overfit and give overly confident predictions on these tasks.The best performance was achieved when the models were trained from a random initialization (contrary to the findings of Wilson et al., 2016a).In free simulation mode RGP performs best of the compared models.This result is expected-RGP was particularly designed to represent and propagate uncertainty through a recurrent process.Our framework focuses on using recurrence to build expressive kernels for regression on sequences.
The suitability of each prediction mode depends on the task at hand.In many applications where the future targets become readily available as the time passes (e.g., power estimation or stock market prediction), the autoregression mode is preferable.We particularly consider autoregressive prediction in the further experiments.

Prediction for smart grid and self-driving car applications
For both smart grid prediction tasks we used LSTM and GP-LSTM models with 48 hour time lags and were predicting the target values one hour ahead.LSTM and GP-LSTM were trained with one or two layers and 32 to 256 hidden units.The best models were selected on 25% of the training data used for validation.For autonomous driving prediction tasks, we used the same architectures but with 128 time steps of lag (1.28 s).All models were regularized with dropout (Srivastava et al., 2014;Gal and Ghahramani, 2016b).On both GEF and self-driving car datasets, we used the scalable version of Gaussian process (MSGP) (Wilson et al., 2015).Given the scale of the data and the challenge of nonlinear optimization of the recurrent models, we initialized the recurrent parts of GP-RNN and GP-LSTM with pre-trained weights of the corresponding neural networks.Fine-tuning of the models was performed with Algorithm 2. The quantitative results are provided in Table 4 and demonstrate that GPs with recurrent kernels attain the state-of-the-art performance.
Additionally, we investigated convergence and regularization properties of LSTM and GP-LSTM models on the GEF-power dataset.The first two charts of Figure 6 demonstrate that GP-based models are less prone to overfitting, even when the data is not enough.The third panel shows that architectures with a particular number of hidden units per layer attain the best performance on the power prediction task.An additional advantage of the GP-layers over the standard recurrent networks is that the best architecture could be identified based on the negative log likelihood of the model as shown on the last chart.Finally, Figure 5 qualitatively demonstrates the difference between the predictions given by LSTM vs. GP-LSTM on point-wise lane estimation (Figure 5a) and the front vehicle tracking (Figure 5b) tasks.We note that GP-LSTM not only provides a more robust fit, but also estimates the uncertainty of its predictions.Such information can be further used in downstream prediction-based decision making, e.g., such as whether a self-driving car should slow down and switch to a more cautious driving style when the uncertainty is high.

Scalability of the model
Following Wilson et al. (2015), we performed a generic scalability analysis of the MSGP-LSTM model on the car sensors data.The LSTM architecture was the same as described in the previous section: it was transforming multi-dimensional sequences of inputs to a two-dimensional representation.We trained the model for 10 epochs on 10%, 20%, 40%, and 80% of the training set with 100, 200, and 400 inducing points per dimension and measured the average training time per epoch and the average prediction time per testing point.The measured time was the total time spent on both LSTM optimization and MSGP computations.The results are presented in Figure 7.
The training time per epoch (one full pass through the entire training data) grows linearly with the number of training examples and depends linearly on the number of inducing points (Figure 7, two left charts).Thus, given a fixed number of inducing points per dimension, the time complexity of MSGP-LSTM learning and inference procedures is linear in the number of training examples.The prediction time per testing data point is virtually constant and does not depend on neither on the number of training points, nor on the number of inducing points (Figure 7, two right charts).

Discussion
We proposed a method for learning kernels with recurrent long short-term memory structure on sequences.Gaussian processes with such kernels, termed the GP-LSTM, have the structure and learning biases of LSTMs, while retaining a probabilistic Bayesian nonparametric representation.The GP-LSTM outperforms a range of alternatives on several sequence-toreals regression tasks.The GP-LSTM also works on data with low and high signal-to-noise ratios, and can be scaled to very large datasets, all with a straightforward, practical, and generally applicable model specification.Moreover, the semi-stochastic scheme proposed in our paper is provably convergent and efficient in practical settings, in conjunction with structure exploiting algebra.In short, the GP-LSTM provides a natural mechanism for Bayesian LSTMs, quantifying predictive uncertainty while harmonizing with the standard deep learning toolbox.Predictive uncertainty is of high value in robotics applications, such as autonomous driving, and could also be applied to other areas such as financial modeling and computational biology.
There are several exciting directions for future research.The GP-LSTM quantifies predictive uncertainty but does not model the propagation of uncertainty in the inputs through a recurrent structure.Treating free simulation as a structured prediction problem and using online corrective algorithms, e.g., DAGGER (Ross et al., 2011), are likely to improve performance of GP-LSTM in the free prediction mode.This approach would not require explicitly modeling and propagating uncertainty through the recurrence and would maintain the high computational efficiency of our method.
Alternatively, it would be exciting to have a probabilistic treatment of all parameters of the GP-LSTM kernel, including all LSTM weights.Such an extension could be combined with stochastic variational inference, to enable both classification and non-Gaussian likelihoods as in Wilson et al. (2016b), but also open the doors to stochastic gradient Hamiltonian Monte Carlo (Chen et al., 2014) (SG-HMC) for efficient inference over kernel parameters.Indeed, SG-HMC has recently been used for efficient inference over network parameters in the Bayesian GAN (Saatchi and Wilson, 2017).A Bayesian approach to marginalizing the weights of the GP-LSTM kernel would also provide a principled probabilistic mechanism for learning model hyperparameters.
One could relax several additional assumptions.We modeled each output dimension with independent GPs that shared a recurrent transformation.To capture the correlations between output dimensions, it would be promising to move to a multi-task formulation.In the future, one could also learn the time horizon in the recurrent transformation, which could lead to major additional performance gains.
Finally, the semi-stochastic learning procedure naturally complements research in asynchronous optimization (e.g., Deisenroth and Ng, 2015).In combination with stochastic variational inference, the semi-stochastic approach could be used for parallel kernel learning, side-stepping the independence assumptions in prior work.We envision that such efforts for Gaussian processes will harmonize with current progress in Bayesian deep learning.eigendecompose a Toeplitz matrix, it can be approximated by a circulant matrix 13 which eigendecomposes by simply applying discrete Fourier transform (DFT) to its first column.Therefore, an approximate eigendecomposition of each T 1 , . . ., T P is computed via the fast Fourier transform (FFT) and requires only O(m log m) time.

A.3 Structure exploiting inference
To perform inference, we need to solve (K SKI + σ 2 I) −1 y; kernel learning requires evaluating log det(K SKI + σ 2 I).The first task can be accomplished by using an iterative scheme-linear conjugate gradients-which depends only on matrix vector multiplications with (K SKI + σ 2 I).The second is done by exploiting the Kronecker and Toeplitz structure of K U,U for computing an approximate eigendecomposition, as described above.

A.4 Fast Test Predictions
To achieve constant time prediction, we approximate the latent mean and variance of f * by applying the same SKI technique.In particular, for a set of n * testing points, X * , we have where KX,X = W K U,U W and KX * ,X = W * K U,U W , and W and W * are n×m and n * ×m sparse interpolation matrices, respectively.Since K U,U W [ KX,X + σ 2 I] −1 y is precomputed at training time, at test time, we only multiply the latter with W * matrix which results which costs O(n * ) operations leading to O(1) operations per test point.Similarly, approximate predictive variance can be also estimated in O(1) operations (Wilson et al., 2015).Note that the fast prediction methodology can be readily applied to any trained Gaussian process model as it is agnostic to the way inference and learning were performed.

Appendix B. Gradients for GPs with recurrent kernels
GPs with deep recurrent kernels are trained by minimizing the negative log marginal likelihood objective function.Below we derive the update rules.
By applying the chain rule, we get the following first order derivatives: The derivative of the log marginal likelihood w.r.t. to the kernel hyperparameters, θ, and the parameters of the recurrent map, W , are generic and take the following form (Rasmussen 13.Wilson et al. (2015) explored 5 different approximation methods known in the numerical analysis literature.
and Williams, 2006, Ch. 5, Eq. 5.9): The derivative ∂K/∂θ is also standard and depends on the form of a particular chosen kernel function, k(•, •).However, computing each part of ∂L/∂W is a bit subtle, and hence we elaborate these derivations below.
Consider the ij-th entry of the kernel matrix, K ij .We can think of K as a matrix-valued function of all the data vectors in d-dimensional transformed space which we denote by H ∈ R N ×d .Then K ij is a scalar-valued function of H and its derivative w.r.t. the l-th parameter of the recurrent map, W l , can be written as follows: Notice that ∂K ij /∂H is a derivative of a scalar w.r.t. to a matrix and hence is a matrix; ∂H/∂W l is a derivative of a matrix w.r.t. to a scalar which is taken element-wise and also gives a matrix.Also notice that K ij is a function of H, but it only depends the i-th and j-th elements for which the kernel is being computed.This means that ∂K ij /∂H will have only non-zero i-th row and j-th column and allows us to re-write (B.7) as follows: Since the kernel function has two arguments, the derivatives must be taken with respect of each of them and evaluated at the corresponding points in the hidden space, h i = φ(x i ) and h j = φ(x j ).When we plug this into (B.6),we arrive at the following expression: (B.9)The same expression can be written in a more compact form using the Einstein notation: where d indexes the dimensions of the h and l indexes the dimensions of W .
In practice, deriving a computationally efficient analytical form of ∂K/∂h might be too complicated for some kernels (e.g., the spectral mixture kernels (Wilson and Adams, 2013)), especially if the grid-based approximations of the kernel are enabled.In such cases, we can simply use a finite difference approximation of this derivative.As we remark in the following section, numerical errors that result from this approximation do not affect convergence of the algorithm.

Appendix C. Convergence results
Convergence results for the semi-stochastic alternating gradient schemes with and without delayed kernel matrix updates are based on (Xu and Yin, 2015).There are a few notable differences between the original setting and the one considered in this paper: 1. Xu and Yin (2015) consider a stochastic program that minimizes the expectation of the objective w.r.t.some distribution underlying the data: where every iteration a new ξ is sampled from the underlying distribution.In our case, the goal is to minimize the negative log marginal likelihood on a particular given dataset.This is equivalent to the original formulation (C.11), but with the expectation taken w.r.t. the empirical distribution that corresponds to the given dataset.
2. The optimization procedure of Xu and Yin (2015) has access to only a single random point generated from the data distribution at each step.Our algorithm requires having access to the entire training data each time the kernel matrix is computed.
3. For a given sample, Xu and Yin (2015) propose to loop over a number of coordinate blocks and apply Gauss-Seidel type gradient updates to each block.Our semi-stochastic scheme has only two parameter blocks, θ and W , where θ is updated deterministically on the entire dataset while W is updated with stochastic gradient on samples from the empirical distribution.
Noting these differences, we first adapt convergence results for the smooth non-convex case (Xu and Yin, 2015, Theorem 2.10) to our scenario, and then consider the variant with delaying kernel matrix updates.

C.1 Semi-stochastic alternating gradient
As shown in Algorithm 1, we alternate between updating θ and W .At step t, we get a mini-batch of size N t , x t ≡ {x i } i∈It , which is just a selection of points from the full set, x.
Define the gradient errors for θ and W at step t as follows: where g T θ and g T W are the true gradients and gT θ and gT W are estimates14 .We first update θ and then W , and hence the expressions for the gradients take the following form: Note that as we have shown in Appendix B, when the kernel matrix is fixed, g W and gW factorize over x and x t , respectively.Further, we denote all the mini-batches sampled before t as Consider the following: a deterministic function of θ t+1 and W t .θ t+1 is being updated deterministically using g t−1 θ , and hence only depends on W t and θ t .Therefore, it is independent of x t , which means that E • Now, consider E[g t W | x [t−1] ]: Noting that the expectation is taken w.r.t. the empirical distribution and K t+1 does not depend on the current mini-batch, we can write: In other words, semi-stochastic gradient descent that alternates between updating θ and W computes unbiased estimates of the gradients on each step.
Remark 1 Note that in case if (∂K t+1 /∂W ) is computed approximately, Lemma 1 still holds since both g t W and E[g t W | x [t−1] ] will contain exactly the same numerical errors.
Assumption 1 For any step t, Lemma 1 and Assumption 1 result into a stronger condition than the original assumption given by Xu and Yin (2015).This is due to the semi-stochastic nature of the algorithm, it simplifies the analysis, though it is not critical.Assumptions 2 and 3 are straightforwardly adapted from the original paper.
Assumption 2 The objective function, L, is lower bounded and its partial derivatives w.r.t.θ and W are uniformly Lipschitz with constant L > 0:

17)
Assumption 3 There exists a constant ρ such that θ t 2 + E W t 2 ≤ ρ 2 for all t.
Lemma 2 Under Assumptions 2 and 3, Proof The inequalities for E θ 2 and E W 2 are trivial, and the ones for E g t Negative log marginal likelihood of a Gaussian process with a structured kernel is a nonconvex function of its arguments.Therefore, we can only show that the algorithm converges to a stationary point, i.e., a point at which the gradient of the objective is zero.
Theorem 1 Let {(θ t , W t )} be a sequence generated from Algorithm 1 with learning rates for θ and W being λ t θ = c t θ α t and λ t W = c t W α t , respectively, where c t θ , c t W , α t are positive scalars, such that (C.18)Then, under Assumptions 1 through 3 and if σ = sup t σ t < ∞, we have Proof The proof is an adaptation of the one given by Xu and Yin (2015) with the following three adjustments: (1) we have only two blocks of coordinates and (2) updates for θ are deterministic which zeros out their variance terms, (3) the stochastic gradients are unbiased.
• From the Lipschitz continuity of ∇ W L and ∇ θ L (Assumption 2), we have: where the last inequality comes from the following (note that g t W := ∇ W L(θ t+1 , W t )): Analogously, we derive a bound on L(θ t+1 , W t ) − L(θ t , W t ): Summing the two obtained inequalities and applying Assumption 1, we get: Therefore, summing the right hand side of the final inequality over t and using (C.18), we have: (C.20) Since the objective function is lower bounded, this effectively means: • Finally, using Lemma 2, our assumptions and Jensen's inequality, it follows that According to Proposition 1.2.4 of (Bertsekas, 1999), we have E g t θ 2 → 0 and E g t W 2 → 0 as t → ∞, and hence where the first term of the last inequality follows from Lemma 2 and Jensen's inequality.

C.2 Semi-stochastic gradient with delayed kernel matrix updates
We show that given a bounded delay on the kernel matrix updates, the algorithm is still convergent.Our analysis is based on computing the change in δ t θ and δ t W and applying the same argument as in Theorem 1.The only difference is that we need to take into account the perturbations of the kernel matrix due to the introduced delays, and hence we have to impose certain assumptions on its spectrum.
Assumption 4 Recurrent transformations, φ W (x), is L-Lipschitz w.r.t.W for all x ∈ X L : Assumption 5 The kernel function, k(•, •), is uniformly G-Lipschitz and its first partial derivatives are uniformly J-Lipschitz: Assumption 6 For any collection of data representations, {h i } N i=1 , the smallest eigenvalue of the corresponding kernel matrix, K, is lower bounded by a positive constant γ > 0.
Note that not only the assumptions are relevant to practice, Assumptions 5 and 6 can be also controlled by choosing the class of kernel functions used in the model.For example, the smallest eigenvalue of the kernel matrix, K, can be controlled by the smoothing properties of the kernel (Rasmussen and Williams, 2006).
Consider a particular stochastic step of Algorithm 2 at time t for a given mini-batch, x t , assuming that the kernel was last updated τ steps ago.The stochastic gradient will take the following form: (C.21)We can define δt W = ĝt W − g t W and uniformly bound δt W − δ t W in order to enable the same argument as in Theorem 1.To do that, we simply need to understand effect of perturbation of the kernel matrix on gt W .
Proof The difference between K −1 t−τ and K −1 t is simply due to that the former has been computed for W t−τ and the latter for W t .To prove the bound, we need multiple steps.
• First, we need to bound the element-wise difference between K t−τ and K t .This is done by using Assumptions 4 and 5 and the triangular inequality: • Next, since each element of the perturbed matrix is bounded by 2GLτ λσ, we can bound its spectral norm as follows: which means that the minimal singular value of the perturbed matrix is at least σ 1 = γ − 2GLτ λσ √ N due to Assumption 6.
• The spectral norm of the expression of interest can be bounded (quite pessimistically!) by summing up together the largest eigenvalues of the matrix inverses: .
Each element of a matrix is bounded by the largest eigenvalue.
Using Lemma 3, it is straightforward to extend Theorem 1 to Algorithm 2.
Theorem 2 Let {(θ t , W t )} be a sequence generated from Algorithm 2 with learning rates for θ and W being λ t θ = c t θ α t and λ t W = c t W α t , respectively, where c t θ , c t W , α t are positive scalars, such that (C.23)Then, under Assumptions 1 through 6 and if σ = sup t σ t < ∞, we have lim t→∞ E ∇L(θ t , W t ) = 0.
• Estimation of the nearest front vehicle position: (x, y) coordinates in the car-centric frame.
-Targets: x and y coordinates.

Appendix E. Neural architectures
Details on the best neural architectures used for each of the datasets are given in Table 5.
Figure 1: (a) Graphical representation of a recurrent model (RNN/LSTM) that maps an input sequence to a target value in one-step-ahead prediction manner.Shaded variables are observable, diamond variables denote deterministic dependence on the inputs.(b) Graphical model for GP-LSTM with a time lag, L = 3, two training time points, t = 3 and t = 4, and a testing time point, t = 5.Latent representations are mapped to the outputs through a Gaussian field (denoted in red) that globally correlates predictions.Dashed variables represent data instances unused at the given time step.(c) Graphical representation of a GP with a kernel structured with a parametric map, φ.

Figure 2 :
Figure 2: Left: Visualization of the GEF-power time series for two zones and the cumulative load with the time resolution of 1 day.Cumulative power load is generally negatively correlated with the temperature measurements on all the zones.Right: Visualization of the GEF-wind time series with the time resolution of 1 day.

Figure 3 &
Figure 3 & Table 2: Left: Train and test routes of the self-driving car in the ENU coordinates with the origin at the starting location.Arrows point in the direction of motion; color encodes the speed.Insets zoom selected regions of the routes.Best viewed in color.Right: Summary of the data collected on the train and test routes.

Figure 4 :
Figure4: Two charts on the left: Convergence of the optimization in terms of RMSE on test and NLML on train.The inset zooms the region of the plot right beneath it using log scale for the vertical axis.full and mini denote full-batch and mini-batch optimization procedures, respectively, while 16 and 64 refer to models with the respective number of units per hidden layer.Two charts on the right: Test RMSE for a given architecture trained with a specified method and/or learning rate.
Point-wise predictions of the lanes made by LSTM (upper) and by GP-LSTM (lower).Dashed lines correspond to the ground truth extracted from video sequences and used for training.LSTM (upper) and by GP-LSTM (lower) position predictions of the lead vehicle.Black markers and dashed lines are the ground truth; blue and red markers with solid lines correspond to predictions.

Figure 5 :
Figure 5: Qualitative comparison of the LSTM and GP-LSTM predictions on self-driving tasks.Predictive uncertainty of the GP-LSTM model is showed by contour plots and error-bars; the latter denote one standard deviation of the predictive distributions.

Figure 6 :Figure 7 :
Figure 6: Left to right: RMSE vs. the number of training points; RMSE vs. the number model parameters per layer; NLML vs. the number model parameters per layer for GP-based models.All metrics are averages over 5 runs with different random initializations, computed on a held-out set.

Table 1 :
Statistics for the data used in experiments.SNR was determined by assuming a certain degree of smoothness of the signal, fitting kernel ridge regression with RBF kernel to predict the targets from the input time series, and regarding the residuals as the noise.Tasks with low average correlation between inputs and targets and lower SNR are harder prediction problems.
out a number of experiments that help to gain empirical insights about the convergence properties of the proposed optimization procedure with delayed kernel updates.Additionally, we analyze the regularization properties of GP-RNN/LSTM and compare them with other techniques, such as dropout.Finally, we apply the model to the problem of lane estimation and lead vehicle position prediction, both critical in autonomous driving applications.

Table 3 :
Mattos et al. (2015)of the models in terms of RMSE on the system identification tasks.The averages were computed over 5 runs; the standard deviation is given in the parenthesis.Results for the RGP model are as reported byMattos et al. (2015), available only for the free simulation.

Table 4 :
Average performance of the best models in terms of RMSE on the GEF and Car tasks.The averages were computed over 5 runs; the standard deviation is given in the parenthesis.

Table 5 :
Summary of the feedforward and recurrent neural architectures and the corresponding hyperparameters used in the experiments.GP-based models used the same architectures as their non-GP counterparts.Activations are given for the hidden units; vanilla neural nets used linear output activations.* rec_dropout denotes the dropout rate of the recurrent weights * Each layer consisted of the same number of units given in the table.*